# R1: 0   -(b)-> T
# R2: T   -(d)-> 0
# R3: T+V -(e)-> Y+V
# R4: Y   -(u)-> p*V
# R5: V   -(c)-> 0

source('subvolume.R')

funk.setup <- function(n, innoc, param_c)
{
	params <- list(
		b=1, d=0.1, e=0.001, u=0.5, p=1000, c=param_c  #c=10
	)
	mV <- ifelse(n > 1, 0.9, 0)

	# transition hazards
	transitions <- function(X) {
		with(c(X, params), {
			c(b, d*X[1], e*X[1]*X[3], u*X[2], c*X[3])
		})
	}

	# stoichiometry matrix
	V <- with(params, {
		matrix(c(
			+1, 0, 0,
			-1, 0, 0,
			-1,+1, 0,
			 0,-1, p,
			 0, 0,-1
		), ncol=3, byrow=TRUE)
	})

	T0 <- params$b / params$d
	X0 <- matrix(c(rep(T0,n*n), rep(0,2*n*n)), nrow=n*n)
	X0[(n+1)*(n+1)/2-n,3] <- innoc

	ctx <- new.env()
	spatial.setup(ctx, X0, transitions, V, n, n, c(0,0,mV))
	ctx$out <- matrix(c(0, as.vector(X0)), 1, 1+3*n*n)
	ctx
}

funk.step <- function(ctx)
{
	oldt <- ctx$t

	spatial.step(ctx)

	if(ctx$t < oldt)  # HACK
		ctx$t <- oldt
	else
		ctx$out <- rbind(ctx$out, matrix(c(ctx$t, as.vector(ctx$X)), 1, 1+3*ctx$W*ctx$H))
}

funk.run <- function(n, innoc, param_c)
{
oldt <- 0
	ctx <- funk.setup(n, innoc, param_c)
	while(ctx$t < 100)
{
		funk.step(ctx)
if(ctx$t > oldt+10)
{
	print(ctx$t)
	oldt <- ctx$t
}
}
	ctx$out
}
